Critical temperature for fermion pairing using lattice field theory 
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Dilute gases of 2-component fermions are of great interest in atomic and nuclear physics. When 
interactions are strong enough so that a bound state is at threshold, universal behavior is expected. 
Lattice field theory provides a first principles approach to the study of strongly interacting systems 
such as this through Monte Carlo simulation. Results of exploratory simulations are presented 
here. In particular, the finite temperature phase transition between superfluid and normal states 
is studied. We present first results for the critical temperature T c and describe the future work 
necessary to determine T c as a function of interaction strength. 

PACS numbers: 03.75.Ss, 11.15.Ha, 34.50.-s, 74.20.-z 
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The ability to trap and cool bosonic atoms so they un- 
dergo Bose-Einstein condensation (BEC) has created a 
unique route toward understanding fundamental proper- 
ties of quantum many-body systems. Cooling fermionic 
atoms to the same temperatures presents additional chal- 
lenges due to Pauli exclusion, but these are being over- 
come 0,0. 

In contrast to BEC, in which a phase transition sepa- 
rates classical from quantum behavior, the quantum na- 
ture of Fermi gases appears gradually as the temperature 
reaches the Fermi temperature Tp. A phase transition is 
expected for fermions which have even the slightest at- 
traction due to the Cooper instability, but this transi- 
tion will occur far below Tp. In the weak coupling limit, 
the system will form a Bardeen, Cooper, and Schrieffer 
(BCS) superfluid below the critical temperature which 
is exponentially small: T c cx Tp exp(— ir/(2kp |a|)) Q, 
where a is the 2-body S wave scattering length The 
Fermi momentum and temperature are defined through 
the number density n: kp = (6n 2 n/ g) 1 ^ 3 and If = kp/2; 
g is the degeneracy. (We will use the natural units 
h = M = Ub = lj where M is the fermion mass.) Weak 
coupling implies kp\a\ <C 1. 

On the other hand, if the attraction between the 
fermions becomes large, a — > — oo and a bound state 
appears at threshold. As 1/a > increases, the fermions 
pair to form bosonic molecules which condense at a tem- 
perature which is experimentally attainable, T c w Q.ZTp. 
It is conjectured that T c will be close to Tp even at 
1/a = ||. Much progress is being made toward produc- 
ing and studying atomic gases in this regime, [E El SHI- 

Furthermore, 1/a — is a special point. In this case 
the only finite physical length scale is the mean inter- 
particle spacing n" 1 / 3 , and the details of the interaction 
between fermions are irrelevant. All physical quantities 
scale like appropriate powers of n times universal con- 
stants, applicable to atomic and nuclear systems. Theo- 
retical study of this limit is difficult because it occurs at 
strong coupling. 

Numerical calculations have been done for this sys- 
tem at zero temperature using fixed-node Greens func- 



tion Monte Carlo and fixed-node diffusion Monte 

Carlo 0] • Within the context of the fixed node approx- 
imation and the assumption that the initial trial wave- 
function has a nonzero overlap with the physical ground 
state, results for the energy-per-particle and pairing gap 
were obtained. These calculations are done at zero tem- 
perature with fixed number of particles, up to 66 so far. 

This work employs a very different Monte Carlo 
method based on lattice field theory [l^ . the same ap- 
proach widely used to compute nonperturbative quanti- 
ties in QCD 0]. It has several advantages over previ- 
ous methods. There is no sign problem in this formu- 
lation, so the fermions can be treated exactly. No as- 
sumptions need to be made regarding the ground state. 
We work in the grand canonical ensemble, so the ther- 
modynamic limit is identically the infinite volume limit. 
Lastly, the theory is naturally formulated for finite tem- 
perature studies. 

Below we present the first numerical investigation of 
the fermion pairing phase transition. The method is de- 
scribed as pedagogically as possible within the length 
constraints. Computations are performed along several 
paths in parameters space, and we observe a sharp de- 
crease in the order parameter for fermion pairing as the 
temperature is increased. These first calculations demon- 
strate the feasibility of using lattice field theory to locate 
and compute the critical temperature for fermion pair- 
ing. Theoretical errors are discussed and can be system- 
atically removed. 

We begin by considering the Hamiltonian for 2 species 
of equal mass particles, Ni and N2 in number: 
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If scattering between particles is limited to low energies, 
then the dominant contribution is the S wave. Further- 
more, if the potential vanishes quickly at long distances 
an effective range expansion can be made for the scatter- 
ing amplitude: (—1/a + \k 2 R + . . . — ik)^ 1 , where a is 
the scattering length and R is the effective range of the 
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potential. In a dilute gas kR « n 1//3 i? < 1, so scatter- 
ing can be completely described by a; the short-distance 
details of the potential are unimportant. The potential 
can be replaced by a local interaction, given in second- 
quantized form by ^Co(ipip) 2 , where tp — (ipi,ip 2 ) T and 
V> = (V'i,V'2) are the fermion fields. Negative Co corre- 
sponds to an attractive interaction. Effective field theo- 
ries have been developed which show how to systemati- 
cally correct the short-distance physics should the need 
arise [HE!- 

We work in the grand canonical ensemble and rewrite 
the partition function Z = Trexp[— f3(H — [iN)] as a path 

integral Z = J VipVip exp (- f£ dx± J d 3 x Cj . The 

fermion fields have anti-periodic boundary conditions in 
the imaginary time (3:4) direction. In order to eliminate 
the Grassmann- valued fields, we must make the Lagrange 
density quadratic in them. To do so an auxiliary scalar 
field (f> is introduced (a la Hubbard-Stratonovich) 
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where to 2 = — Cq . 

The Lagrangian is invariant under a global U(l) trans- 
formation ip — > exp(ia)ip, ip — y ipexp(-ia). If we want 
to study the spontaneous breaking of this symmetry in 
a finite volume, we must induce an explicit breaking and 
compute observables in the limit of its removal. Thus 
we add to the Lagrangian a term proportional to a new 
parameter J: 



C[J] 
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£ + -(J^ T a 2 ip + J*il>a 2 iP T ) . (3) 



Without loss of generality we take J to be real and de- 
fine the pairing condensate £ = (ip 1 a 2 i> + ^p<y 2 ip T ) /2. 
limj^oOimy^oo E) is an order parameter for the break- 
ing of the U(I) symmetry and the emergence of a nonzero 
pairing condensate. For the sign convention of J > 
will induce £ < 0. 

C is discretized on a lattice with V spatial sites and 
N t sites in imaginary time; we denote the spatial and 
temporal lattice spacings by b s and b t . To perform simu- 
lations, path integrals must be manipulated into the fol- 
lowing form |f2T | : Combining ip and ■;/> into a 4-component 
fermion, ^ T = (tp T \\p{ia2) T ) 1 the action can be written 
as ^> T A^, where the WN t x WN t matrix A is 
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and is the only 2VN t x 2VN t block to contain nontrivial 
spacetime terms. The only coupling between species is 
through the ia 2 blocks. We also introduce the anisotropy 
£ = b 2 s /b t . The result of the integration is 

Z= I V<j> det \\J\ 2 +K J{ K] e- s *- = [ V<f)W{<f>} (5) 



where K = K®I 2 defines the VN t x VN t matrix K. The 
important point is that W, the integrand of Z, is strictly 
nonnegative, so it can be interpreted as a probability 
weight corresponding to a particular field configuration. 

The generating functional for fcrmionic observables is 
obtained by adding a fermionic source term 3^ to C and 
integrating over 'J. The result is 



m = J 



V<pW[4>] exp 
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The expectation value of a general fermion bilinear is 
generated by derivatives of Z[3] as follows: 



-Tr BA~ 



(7) 
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The interesting bilinears in this work are the fermion den- 
sity n and the pairing condensate E, for which 



i<J 2 S^ 
ia 2 S_ 
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where S± denotes a temporal shift, S±ip x = ipx±e4,- In 
practice a number of random sources N src are used to 
evaluate the trace. Let r)i(x, a) be a Gaussian-distributed 
random complex number at site x; a labels the 4 in- 
ternal \1/ indices, and i runs from 1 to N src . Then 
TiBA^ 1 = J2i 1 l\BA^ 1 ri l /N src . The fermion density can 
also be computed through n = m 2 (</>} = m 2 J Vcf) <j) W [4>] ; 
this provides a check of the code used to compute J7J). 

As stated earlier, this lattice field theory of a homo- 
geneous, dilute gas of 2-component fermions starts from 
first principles, and uncertainties can be removed sys- 
tematically. The thermodynamic limit is reached sim- 
ply by increasing the volume of the system. The con- 
tinuum limit is more nuanced. To remove the lattice 
spacing b s while maintaining constant physics, one keeps 
n 1 / 3 a held fixed; then all physical length scales are re- 
lated to rt~ 1//3 . The continuum limit is the limit where 
n^ 1 ' 3 3> b s , achieved by tuning lattice parameters so 
that (b s /a) -> and /i6 2 /£ -> 0. 

There are 3 physical quantities we wish to control: n, a, 
and T. We do so by tuning fi, to 2 , and the size of the 
lattice in imaginary time, varying either iV t or £. A cal- 
culation in effective field theory matches the scattering 
length to the lattice parameters m 2 and £, modulo finite 
volume and spacing corrections |f2T |: 
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pi = 2sin(pi/2) and the BZ integration/summation is 
over available momenta in the Brillouin zone. The tem- 
perature is tuned by varying £ 



T 



2/3 iV t (n & 3 ) 2 / 3 



(10) 
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where no denotes the zero temperature density. One 
expects the chemical potential /j to strongly control n. 
However, the exact mapping of lattice parameters to 
physical quantities is not is not known a priori and must 
be determined via simulation before any precise physics 
results can be obtained. This paper is the first voyage 
into this parameter space. 

The equivalence of the quantum statistical path inte- 
gral JSJl to a classical statistical partition function sug- 
gests that numerical methods for simulating spin systems 
like the Ising model can be used to study this dilute 
Fermi gas. Additionally we can draw from the exper- 
tise developed for lattice QCD. For any configuration of 
field values C = {4> x }i the contribution of any observable 
Q[4>] is weighted by the nonnegative factor W[<j>]. One 
uses importance sampling to generate an ensemble of 4 
dimensional field configurations which give the greatest 
contributions to (Q) = Z^ 1 fV(j)Q[(j)] W[4>]. Having ob- 
tained a sample of N configurations which are distributed 
according to W[t/>], observables are estimated by their en- 
semble average Q = Q[C]/iV. 

In this work, the hybrid Monte Carlo (HMC) algorithm 
(lif is used to generate samples of scalar field configura- 
tions. One can think of the sequence of configurations as 
successive snapshots of scalar fields evolving in an artifi- 
cial time r, referred to as Monte Carlo (MC) time (dis- 
tinct from Euclidean time X4). Given some initial config- 
uration Co = C(to), a trial configuration C (to + At) is 
obtained by evolving the 4D scalar fields forward in r us- 
ing classical equations of motion. The process Co — > C 
is called a trajectory and is done in discrete steps. The 
change in the HMC Hamiltonian AH is computed, and 
the new configuration is accepted as C\ with probabil- 
ity min(l, e An ). The effects of the fermions are included 
by rewriting the determinant as a path integral over a 
complex scalar field x with action x*{\J\ 2 + K^K)~ 1 \- 
The x fi e ld is held fixed during the classical evolution. 
At the end of a trajectory it is updated by a heatbath 
step: x = iJ*Vi + K^r/2, where 771, 772 are Gaussian noise 
vectors. 

In some regions of parameter space, it is possible that 
the attraction between fermions is so strong that the lat- 
tice is densely packed with a fermion at every site, even at 
/i = 12]. This lattice phase would be separated from 
the dilute, continuum-like phase by a first order tran- 
sition, preventing extrapolation to zero lattice spacing. 
Happily, it is straightforward to check that nb 3 <C 1 for 
/i — in the interesting region of (m 2 ,£) space, namely 
where b s /a = 0|- For example, on a lattice with 
V = 6 3 and N t = 12, nb 3 s = 0.0073(6) when fi = at 
m 2 = 0.1456, £ = 1 (corresponding to 1/a = in infinite 
volume). Due to finite volume effects n > 0. The lesson is 
that these lattice parameters correspond to simulation of 
a dilute continuum- like Fermi gas, not a densely packed 
lattice of fermions. Consequently, we can study thermo- 
dynamics at ijl > confident that the continuum limit 
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FIG. 1: The fermion condensate £ (left) and the density 
n in lattice units (right) as functions of external source for 
m 2 = 0.175. Different symbols indicate different values of £, 
corresponding to different temperatures. Points at J > are 
simulation data, and those at J = are the results of linear 
extrapolations (solid lines). Some data points are displaced 
horizontally for clarity. 



can be reached in the manner described above. 

The results presented below were obtained on a V = 8 3 
lattice with fj,bt — 0.4. With Nt — 16, m 2 was set to 1 of 
4 values, and then 4-5 values of £ were used to locate the 
phase transition. At each (m 2 ,£) pair, separate simula- 
tions were performed setting J = 0.07, 0.1, and 0.14 in 
order to extrapolate to J — 0. Each simulation was run 
for 4000 HMC trajectories with 32, 25, and 20 steps per 
trajectory for respective values of J. Observables were 
computed (with N SIC — 10) at 20 trajectory intervals for 
a total of 200 measurements per simulation. Dropping 
the first 50 measurements was sufficient to ensure the 
sample had equilibrated. To estimate correlations be- 
tween successive measurements, the data were binned in 
groups of 2, 5, 10, and 20; iVbin = 5 sufficed to account 
for correlations. The statistical errors are certainly much 
smaller than the finite volume and lattice spacing uncer- 
tainties of this exploratory study, so we do not discuss 
them further. 

Fig. ^ shows the J — > extrapolation of E for m 2 = 
0.175; the extrapolations for the other values of m? are 
similar. A more accurate extrapolation, exploiting an 
effective field theory analysis can be made once more data 
are obtained; however, in this work we emphasize the 
sharp decrease in E| j^q. Corrections to linearity will not 
qualitatively change our exploratory conclusions about 
locating T c . 

In order to convert the observation of a phase transi- 
tion at a critical £ into a critical temperature using l|10|l . 
we need the zero temperature fermion density tiq. In 
this work we assume that T c /4 is sufficiently low com- 
pared to T c that we can compute no with Nt = 48 (and 
J = 0.14). We find n Q b 3 s = 0.2612(6) and 0.2893(7) for 
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FIG. 2: The fermion condensate E, after extrapolating 
J — > 0. As £ increases, so does T, and the system goes from 
superfiuid to normal. Lines connect the points merely to guide 
the eye. 
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FIG. 3: Critical temperature vs. inverse scattering length. 
The solid line shows T c in the BCS regime and the bold bro- 
ken line shows the T c for BEC of di-fermion molecules. In 
the center is the present, exploratory result for T c /Tf- The 
horizontal error bar reflects the sensitivity of a on £ over the 
critical region. 



(m 2 ,£) = (0.155,1.2) and (0.175,1.3), respectively. (The 
quoted uncertainties represent the statistical error only.) 
The densities at T c /4 are close to those at T c , e.g. the 
corresponding nb 3 s = 0.2577(12) and 0.2877(7). 

Little dependence of n on N t or J is apparent in the 
finite temperature data. Results for n with Nt = 16 and 
m? = 0.175 are shown in Fig. ^ plots for other values of 
m 2 are similar. 

Using these low temperature densities in (|10|) leads to 
T C /T F = 0.035 ± 0.004 and 0.036 ± 0.004 for m 2 = 0.155 
and 0.175. It turns out that b s /a changes significantly 
as £ is tuned through the transition region. With m 2 = 
0.155, b s /a goes from -0.10 at £ = 1.1 to 0.10 at £ = 1.2, 



and with m 2 = 0.175, b s /a varies from —0.11 at £ = 1.2 
to 0.08 at £ = 1.3. In order to determine T c precisely for 
fixed values of 1/a, simulations should hold £ fixed and 
vary /x across the transition. This work is in progress. 

Fig. compares the results of this work with calcula- 
tions of T c in the BCS 0] and BEC regimes. Since the 
different curves in Fig. [21 correspond to the phase transi- 
tion at slightly different lattice spacings (i.e. slightly dif- 

1 /3 

ferent n b s ) but over the same large ranges of l/(n a), 
only one data point is plotted. Reduction of lattice spac- 
ing and volume errors will allow T c to be studied more 

1/3 

precisely in the region — 1 < n a < 1. 

In conclusion, these initial results demonstrate both 
the advantages and the challenges of using lattice field 
theory to study fermion pairing at strong coupling. The 
critical temperature is straightforward to find; we showed 
the condensate vanishing across the finite T transition for 
4 paths through parameter space. The scattering length 
varies quickly along those paths, so different directions in 
parameter space must be used to compute T c while hold- 
ing a fixed. The effects of lattice volume and spacing 
need to be studying systematically and will require sig- 
nificant computational resources. This exploratory work 
represents the first leg of the journey to obtain a first 
principles calculation of T c in the universal regime. 
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